#!/usr/bin/perl -w
use strict;
use DBI;
use DBD::Pg;
use Getopt::Long;
use Digest::MD4;
use FindBin;
use lib ("$FindBin::Bin/..");
use parse_swissprot;

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#	Usage

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
my $usage = <<'USAGE';

USAGE:

    panda_pubmed_from_sprot.pl  --dbname STRING [--verbose] [FILE]

    Mandatory arguments:

         --dbname STRING    the postgreSQL database
    
    Optional arguments:
        
        [--verbose]         dots printed to STDERR every 2000 records
        [--user STRING]     the postgreSQL user name
        [--host STRING]     the postgreSQL database server
        [--password STRING] the postgreSQL database password
        [--error_file FILE] holds parsing errors

        [FILE]              The SwissProt file you wish to process e.g. swissprot.dat
                            This can be specified on the command line or via STDIN 
                            

        [--help]
                            
RESULTS:
    

    This script extracts pubmed IDs from SwissProt and inserts them into the Panda
	 database

USAGE


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#    Command line options

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
# options variables
my $user			= 'postgres';
my $password		= '';
my $host			= '';

my $help			= '';
my $dbname			= '';
my $error_file		= '';
my $verbose;


#output
GetOptions(	'verbose'		=> \$verbose,
			'user=s'		=> \$user,
			'password'		=> \$password,
			'host=s'		=> \$host,
			'dbname=s'		=> \$dbname,
			'error_file=s'  => \$error_file,
			'help'	   		=> \$help);

die $usage if ($help);
die "\nERROR:\n\tMissing Arguments [dbname]:\n" .$usage
				unless ($dbname);

my $md4 = new Digest::MD4;

	

	

	

	
#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


#	uniquify

#			removes duplicate entries in an array

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

sub uniquify(\@)
{
	my %hash;
	@hash{@{$_[0]}} = ();
	@{$_[0]} = keys %hash;	
	return @{$_[0]};
}	






#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#	per_record_call_back

#		parse_pfam::parse_records will call this routine once per PFAM record

#		Description
#			A) Process SQ fields for the protein sequence

#			B) Get Pubmed IDs associated with this sequence (via panda numbers).


#		Paramters:
#			1) record data: a reference to a hash whose keys are field names
#							 e.g. $pfam_data->{'DE'} is the description field
#							 these are all arrays (ie. $pfam_data->{'DE'}[0..n])
#							 except for $pfam_data->{'SQ'} which contains the sequence
#			2) file position of the record start
#			3) record number (indexed from 1) of the current record
#			4) user defined data. In this case, it is an array of
#								taxid_to_name lookups
#								file handles for pubmed, omim codes, variants

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
sub per_record_call_back($$$$)
{
	#$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
	#	FASTA Sequence

	my ($pfam_data, $recstart, $reccnt, $custom_data) = @_;
	my ($verbose			,
		$cnt_pubmed,
		$run_cmd			) = @$custom_data;

	#print progress dots
	if ($verbose && $reccnt % 2000 == 0)
	{
		print STDERR "\n\t" if (($reccnt / 2000) % 50 == 0);
		print STDERR '.';
	}

	# make sure all required fields are present
	for (qw (SQ))
	{
		die "\nFatal Error:\n".
			"The field '$_' is missing ".
			"in record # $reccnt beginning on line $recstart\n"
			unless exists $pfam_data->{$_};
	}
	
	return unless $pfam_data->{'RX'};
	
	#-------------------------------------------------------------------------------------
	# get pi numbers
	my @digest = unpack("L*", $md4->hash($pfam_data->{'SQ'}));
	#
	#-------------------------------------------------------------------------------------
	
	#-------------------------------------------------------------------------------------
	# pubmed ID
	# RX
	# save unique pubmed IDs
	my @PubMed;
	push(@PubMed, $1) if(/PubMed=(\d+)/) foreach (@{$pfam_data->{'RX'}});
	uniquify(@PubMed);
	
	$$cnt_pubmed += @PubMed;
	
	# for each pubmed ID, print out 128-bit hash of the sequence to which it is linked
    print $run_cmd join ("\t", @digest, "$_\n") foreach (@PubMed);
	#-------------------------------------------------------------------------------------

}


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#	main

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

print STDERR "\tParsing Swiss Prot file...\n\t[" if ($verbose);

my $fh_SwissProt = *STDIN;
if (@ARGV)

my $cnt_pubmed		= 0;
my $reccnt          = 0;

			'user=s'		=> \$user,
			'password'		=> \$password,
			'host=s'		=> \$host,
			'dbname=s'		=> \$dbname,

my $exec = "| $FindBin::Bin/../../add_sequence_linked_data ";
$exec .= "--dbname $dbname ";
$exec .= "--user $user " if length($user);
$exec .= "--host $host " if length($host);
$exec .= "--password $password " if length(password);
$exec .= "--error_file $error_file " if length(error_file);
$exec .= "--verbose " if verbose;


my $first_pass = true;
for my $file_name(@ARGV)
{
	my $run_cmd = $exec;
	$run_cmd .= "--remove_first " if $first_pass;
	open RUN_CMD, $run_cmd 
				or die "Error\n\tCould not find or run add_sequence_linked_data\n$@\n";
	my $cnt_pubmed_in_file		= 0;
	open SP_FILE, $file_name or die "Error\n\tCould not open $file_name\n$@\n";
	$reccnt += parse_records(*SPFILE, 					# SwissProt file handle
							 @{[qw (SQ RX)]},			# desired sprot fields
							 \&per_record_call_back,	# call back function
							 [$verbose, \$cnt_pubmed_in_file, *RUN_CMD] );
	$cnt_pubmed += $cnt_pubmed_in_file;
	close RUNC_CMD;
}

# swiss prot seems to use taxids which are no longer present in NCBI!
print STDERR "]\n\t$reccnt\tSwissProt records were parsed.\n" if ($verbose);

# print summary
if ($verbose)
{
	print STDERR "\t", $cnt_pubmed, "\tpubmed IDs parsed.\n";
	
}
print STDERR "\n";


